function [LL, varargout] = loglikelihood_by_blocks(modelId, data, params, takeNegative)
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% This function computes the same thing as loglikelihood but avoids large RAM requirements by computing things by block, running
	% a for loop across blocks to sum up results.
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Inputs:
	% modelId:			       integer (1 => BinomialLogistic, 2 => Poisson)
	% data:			           object:
	% 	.LL_offset:                scalar
	%   .blocks:                   cell(NumBlocks, 1)
	% params:                  NumX x 1
	% takeNegative:            boolean
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Outputs:
	% LL:            		scalar
	% grad:					NumParams x 1
	% hessian:				NumParams x NumParams
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	NumBlocks = length(data.blocks);

	%%% Call function loglikelihood for each block + sum up results to grand total each time
	if nargout == 1
		LL = 0;
		for bb = 1:NumBlocks
			LL = LL + loglikelihood_memory_efficient(modelId, data.blocks{bb}, params, takeNegative);
		end
	end
	if nargout == 2
		LL = 0; grad = 0;
		for bb = 1:NumBlocks
			[LL_bb, grad_bb] = loglikelihood_memory_efficient(modelId, data.blocks{bb}, params, takeNegative);
			LL   = LL   + LL_bb;
			grad = grad + grad_bb;
		end
		varargout{1} = grad;
	end
	if nargout == 3
		LL = 0; grad = 0; hessian = 0;
		for bb = 1:NumBlocks
			[LL_bb, grad_bb, hessian_bb] = loglikelihood_memory_efficient(modelId, data.blocks{bb}, params, takeNegative);
			LL      = LL      + LL_bb;
			grad    = grad    + grad_bb;
			hessian = hessian + hessian_bb;
		end
		varargout{1} = grad;
		varargout{2} = hessian;
	end
	
	% Add LL_offset at the end
	if ~takeNegative
		LL = LL + data.LL_offset;
	else
		LL = LL - data.LL_offset;
	end
	toc
end
